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Abstract 



•S 

The reduced basis method is a model reduction technique yielding substantial savings 
of computational time when a solution to a parametrized equation has to be computed 
for many values of the parameter. Certification of the approximation is possible by means 
of an a posteriori error estimator, delivering an upper bound on the error. Under ap- 
propriate assumptions, this error bound is computed with an algorithm of complexity 
^ independent of the size of the full problem. In practice, the evaluation of the error esti- 

Q\ mator can become very sensitive to round-off errors. We propose herein an explanation of 

this fact. A first remedy has been proposed in [F. Casenave, C. R. Acad. Sci.]. Herein, 
{vq we improve this remedy by proposing a new approximation of the error estimator using 

t-H the Empirical Interpolation Method. This method achieves higher levels of accuracy and 

requires potentially less precomputations than the usual formula. A stabilized version is 
also derived. The method is illustrated on a simple one-dimensional diffusion problem and 
a three-dimensional acoustic scattering problem solved by a boundary element method. 

1991 Mathematics Subject Classification. 65N15, 65D05, 68W25, 76Q05. 
Keywords. Reduced basis method, a posteriori error estimator, round-off errors, boundary 
element method, empirical interpolation method, acoustics 



Introduction 

In many problems, such as optimization, uncertainty propagation or real-time simulation, one 
has to evaluate an objective function for a large number of values of some parameters. Evaluat- 
ing this objective function often implies solving a parametrized partial differential equation for 
the given parameter value. In an industrial context, one evaluation of the objective function 
can already be a challenging numerical problem. To keep reasonable computational costs, 
various model reduction techniques have been developed to speed up computations, such as 
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Guyan condensation [22] and dynamic condensation [33] in vibration of structures, proper 
orthogonalized decomposition [12], proper generalized decomposition [15], sparse polynomial 
approximation [37], projection-based reduced order models [I], tensor product approximations 
with greedy techniques [8, 31 J . . . We focus in this document on another procedure, the Reduced 
Basis (RB) method [28, 35J. This method has been applied to many kinds of problems, includ- 
ing non-linear problems such as the viscous Burgers equation [40J or the steady incompressible 
Navier-Stokes equations [59] . 

As described in Section [T| the RB method consists in replacing the sequence V B \x h- > 

|_[ 

u^i | — > Q(Ufj,) by the sequence V B /i H > \- > Q(u^). Here, V denotes the parameter set, 

: \i h-> u M the model problem, : fi \-t its lower-dimension approximation, Q(it M ) 
the quantity of interest, and Q(u M ) its RB approximation. More specifically, the RB method 
consists in two steps: (i) A so-called offline stage, where solutions of E^ for well-chosen 
values of the parameter \i are computed. During this stage, N problems of size N are solved 
(with N <C A), and some quantities related to the N solutions are stored; (ii) A so-called 
online stage, where the precomputed quantities are used to solve E^ for many values of ji. 
In this stage, an a posteriori error estimator £(fi) is computed to check the quality of the 
approximation on the quantity of interest. An important feature in the reduced basis method 
is the use of an efficient error estimator, which is an upper bound of the error made while 
considering the reduced solution. This is the so called certification of the error. The notion 



of efficiency is defined in Section 1.4 The error estimator must be as sharp as possible to 
faithfully represent the error. However, as noticed for example in |33|, pp. 148-149], the error 
estimator is subject to round-off errors, especially in compliant academic examples. The 
purpose of this work is an explanation of this fact and the derivation of a new method to 
compute the estimator in an accurate and efficient way using potentially less precomputed 
quantities than the classical formula. 

In Section [TJ we recall the main ingredients of the RB method, namely (i) the construction 
of the reduced problem, (ii) the a posterior error estimator, (iii) the notion of efficiency, (iv) 
the offline stage during which the vectors of the reduced basis are constructed. We then 
explain in Section [2] why the classical formula for computing the a posteriori error estimator 
is ill-conditioned in regard of round-off errors. In Section [3j we present our new procedure 
based on the Empirical Interpolation Method. A stabilized version is also derived, and the 
various procedure to compute the error estimator are compared on a simple one-dimensional 
diffusion problem. In Section [4] we apply this new procedure to a three-dimensional acoustic 
scattering problem. 



1 The reduced basis method 
1.1 The model problem 

We suppose that the problem of interest has the following discrete variational form, depending 
on a parameter [i in a parameter set V- for a finite-dimensional space V of dimension N (with 
JV» 1 resulting, e.g., from discretization), find G V such that 

En : a^iu^ v) = b(v), VdG V, (1) 

where is an inf-sup stable bounded sesquilinear form on V x V, antilinear with respect to 
the right variable, and b is a continuous linear form on V. We work in complex vector spaces 
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in view of our application to acoustic scattering. In what follows, the complex conjugate of 
z G C is denoted z* . We define the Riesz isomorphism J from V' to V such that V/ G V', 
Vu G V, (Jl,u) v = l{u), where (-,-)v denotes the inner product of V with associated norm 

II • ||v- We denote j3u := inf sup } ^- ' }} > the inf-sup constant of % and /L a computable 

ueV ve v IMIvlMlv 

positive lower bound of (5 p. For simplicity, we consider that the linear form b is independent of 
the parameter [i. The extension to /^-dependent b is straightforward. We refer to the discrete 
solution Up as the "truth solution" . We denote Ap the stiffness matrix of size N related to 
an, once a basis of V has been selected. In this basis, we denote B the source term and Up 
the solution, which are vectors of size N. The problem En then consists in solving the linear 
system A lx U tl = B. 



1.2 The reduced problem 

Suppose that a reduced basis, consisting of N solutions Up i of E^, i G {1, N}, has already 
been constructed. To alleviate the notation, we denote Ui the function u^. We will see in 



Section [L5| how the parameters m are chosen. Given a parameter value /i G V, the reduced 
problem is then a Galerkin procedure written on the linear space V = Spanjui, ...,u^} C V: 
find Up G V such that 

Ep:ap(up, Ul ) = b{ Uj ), Vj G {1, ...,N}. (2) 
The approximate solution on the reduced basis is written as 

N 



J^7i(M)«i- (3) 



% = 

In matrix notation, Up = (7i(M))i<»<JV ^ s som ti° n t° * ne linear problem 

= ^- ( 4 ) 

where (^)ij = dp{ui,Uj) and = for 1 < i,j < JV. 

Recalling the exact and approximate quantities of interest Q(up) and Q(u,p), respectively, 
the quality of the approximation for a given \i G V is quantified by the error measure 
— 0C*v)||- Two main cases are generally considered: (i) the so-called general-purpose 
(or compliant) case, where we are interested in the whole solution: Q = Q = Id and || • || = || • ||y ; 
(ii) the so-called goal-oriented (or non-compliant) case, where Q is a linear form on V and 
|| • || = | • |. When we obtain a satisfying approximation with N <C N, the RB strategy is 
successful. The operator Q is consistently built so that ||Q(u M ) — Q{up)\\ vanishes for u = Ui, 

{1 V}. 



1.3 A posteriori error estimator 

In the standard RB method, the a posteriori error estimator is a residual based estimator. In 
what follows, we refer to it simply as error estimator. Since this error estimator is an upper 
bound, it provides a way to certify the approximation made by the reduced basis. 

Property 1.3.1 (General-purpose case). The following error bound holds: For all \i G V , 

\\up - Up\\ < £i((i), (5) 
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where 

:= ^IIG^Hv, (6) 
with G^ the linear map from V to V such that V9mi4 G^u := J (a M (u, •) — &)£ V. 

Proo/. See [331 Section 4.3.2]. □ 

a first 



1.3.1 



In the goal-oriented case, Q £ V' so that choosing Q = Q and using Property 
error estimator is 

Q( Ufl ) - Q(6 M ) < £ go (/i) := ^[iQllvllG^llv (7) 

It is possible to obtain a more accurate approximation of Q{u^) by slightly modifying the 
quantity of interest. We introduce the following dual problem: Find ^GV such that 

E*:a li (w,v fl ) = Q(w), Vu> E V. (8) 

We wrote the dual problem on the same discrete space V, but another space can be considered. 
A reduced basis procedure is also carried out for the problem E$, resulting in an approximation 
Vfj, of V/j,. The modified quantity of interest is then defined as Q(u^) := Q(u^) — (G^u M , v^v, 
where the second term is the so-called dual-based correction of residual estimators. 

Property 1.3.2 (Goal-oriented case). The following error bound holds: For all fj, £ V, 

Q(u) - Q«>| < := IIG^IIvllG^Hv, (9) 

where G^ is the linear map from V to V such that V 3 v \— > G^u := J(a /Jt (-,v) — Q) E V 

and M is a computable lower bound of j3f. = inf sup ] ^ ' )} . Obviously, /3f. = P u if a u is 

* M uev v( z V \\u\\v\\v\\v 

hermitian. 

Proof. See [3 Proposition 23] and [TTJ Proposition 3.1]. □ 

The error estimator in ^ involves the product of two residuals, whereas the error estima- 
tor in ([7]) involves only one residual. Therefore, the order of convergence with respect to the 
residuals is twice for Q as it is for 0. 

In what follows, we mainly focus on the general-purpose case. Extensions to the goal- 
oriented case are straightforward. 

1.4 Efficiency of the RB method 

The notion of efficiency is central to the RB method. 

Definition 1.4.1. The RB method is efficient if (i) the reduced problems are constructed in 
complexity independent of N , and (ii) the error estimator is computed in complexity indepen- 
dent of N . 

The following assumption is sufficient for the RB procedure to be efficient. 
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Definition 1.4.2. The sesquilinear form depends on fi in an affine way if there exist d 
functions : V — > C and d ^-independent sesquilinear forms ai bounded on V x V such 

that 



a^(u, v) = } J a k (n)a k (u, v), Vu, v G V. 



(10) 



fc=i 



In what follows, we always assume that the affine decomposition (10) holds. 



Property 1.4.1. If depends on ji in an affine way, then the RB procedure is efficient. 

Proof, (i) The reduced matrix writes = a k (fi)A k , where (A k )ij := a k (ui,Uj). There- 

fore, provided the d matrices A k and the vector B are precomputed during the offline stage, 
the reduced problems are constructed in complexity independent of N. 

(ii) The operator G M inherits the affine dependence of a M on /x since, Vit G V, 

d d 
G^u = -Jb + a k (fi)Ja k (u, •) = G 00 + a k (fj,)G k u, 



where Goo := 
decomposition, we infer 



k=i 

Jb £ V and G^ri := 



Ja k (u, ■ 



k=l 

G V for all fc G {1, 



N 



G, 



00 



+ ^ a k (n)^i(^)G k Ui 
i=i fc=i 



(11) 

, d} . Using this affine 
(12) 



The scalar product on which the norm in (12) hinges can be expanded to provide another 
formula for the error estimator (see |33[ eq.(4.61)]): 

/ N d 

£ 2 {n) = ft' 1 (G o,Goo)v + 2Re^2^2ji([i)a k (n)(G k Ui,G o) v 



1=1 k=l 



N d 

■ ^2 ^2 7i(^)afc(^)7j(M)"?(A t )(G'fcUi,G/'u j )i 

i,j=l k,l=l 



(13) 



which is computed in complexity independent of N provided that (Goo,Goo)vj (G k U{, Goo)v 
and (G k Ui, G;Uj)v are precomputed during the offline stage, and provided that a lower bound 
fin of the stability constant of is also computed in complexity independent of N (which is 
possible by the successive constraint method, see (2BJ). □ 



Furthermore, an important observation made in [9\, and that will be usefull below, is that 



the formula ( 13 ) defining £ 2 can be rewritten in an equivalent way as 

£ 2 (n) := fin 1 {S 2 + 2Re(s t Xn) + x^Sx^) 1 , 



(14) 



where 5 := ||Goo||v; s an d x^ are vectors in C with components sj := {G k Ui,Goo)v and 
{xn)i '■= a k{v)li{v) ■> an d S is a matrix in <C dN > dN with coefficients Sj y j := (G k Ui,GiUj)v 
(with / and J re-indexing respectively (k,i) and 1 < k, I < d, 1 < i, j < N). The t 

superscript denotes the transposition. The vector s and the matrix S depend on the reduced 
basis functions {^Ik^^t but are independent of and the vector x M depends on the RB 
approximation via the coefficients 7i(/i). Notice that the term between parenthesis in the 
right-hand side of (14) is a multivariate polynomial in x„ of total degree 2. 
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1.5 The offline stage 

Fix a discrete subset of parameters Ptriai C V . In the offline stage, the parameters \ii (from 
which the reduced basis is constructed using the solution to E^) are chosen by a greedy 
algorithm as elements of Ptriai- We denote Select the set of these selected parameters. Let 
tol be the value of the desired certification for the error. The greedy algorithm proceeds 
as described in Algorithm [T] (see [331 Section 3.3]). At each step of the algorithm, the new 
quantities a^{ui,Uj) and b(v,j) are computed and stored, as well as the new components of the 



vector s and of the matrix S to be used in the formula £2 ( 14 ). This task, as that of evaluating 
GoO) typically requires inverting the stiffness matrix in V by solving, for all k £ {1, ...,d} and 
all i £ {1, N}, the variational problem: find Wik £ V such that 

Ea,k ■ (wi,k,v)v = a k (ui,v), \/v £ V. (15) 

Then, G^Ui = Wi : k can be computed. The computation of (GkUi, GiUj)v follows from the 
solutions of Ea t k an d Egjj- Since the error estimators are evaluated using the formula 
£2(^)1 fJ- £ Atrial, with the current state of the reduced basis, finding the maximum of the error 
estimator on Ptriai is of complexity independent of N. This allows one to consider very large 
sets Ptriai without increasing too much the complexity of the whole offline procedure. Finally, 
notice that for goal-oriented RB, each of the steps in Algorithm [T] has to be carried out for 
the adjoint problem ^ as well. 

2 Round-off errors and online certification 



In this section, we explain why the online error estimator (14) may be sensitive to round-off 
errors. 

2.1 Elements of floating point arithmetic 

In a computer, real numbers are represented by a finite number of bits, called floating-point 
representation. Current architectures are optimized for a format used by a large majority of 
softwares: IEEE 754 double-precision binary floating-point format. Elements of floating-point 
arithmetic can be found in [J, Chapter 1] and in [20] . A real number is represented in scientific 
notation as 

(-l) sign mantissa x 2 ex P° nent - 1023 . 

This format contains 64 bits of data: 1 bit for the sign, 52 bits for the mantissa, and 11 bits 
for the exponent. In base 10, real numbers are represented with roughly 16 digits, and the 
minimum and maximum exponents are —308 and 308. Let x be a given real number. We 
denote fl(x) the floating point representation of x. In what follows, we only consider reals 
such that 10~ 308 < \x\ < 10 308 . 

Proposition 2.1.1. When a real number is rounded to the closest floating-point number, the 
relative error on its floating-point representation is bounded by a number, e, called the machine 
precision: 

y x G [_io 3 ° 8 , -IO- 308 ] U [IO- 308 , 10 308 ], 



fl(x) 



< e. 



In base 10, e = 5 x 10~ 16 (see [201 Section 1.2]). 

Consider an operation (g) (that can be a summation, a subtraction, ...). To treat the 
operation (g> between two reals x and y, the computer returns fl{fl{x)®fl{y)). In some cases, 
this quantity can be very different from the exact result x (8> y. 
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19. 
20. 



23. 
24. 

25. 
26. 
27. 



30. 
31. 
32. 
33. 
34. 
35. 
36. 



{Mi} 



[Initialization of Select] 
[Constant term of £2] 



Algorithm 1 Offline stage of the RB method 

1. Set N = 1 

2. Choose m € Ptrial randomly and set Select 

3. Compute and store (Goo,Goo)v 

4. Compute u\ solution of E^ 

5. Set Vi = Spanjui} [Initialization of the reduced basis] 

6. Compute and store b (u±) [First coefficient of the reduced right-hand side B] 

7. for all k in {1, d} do 

8. Compute G^ui as the solution of Eqi^ [Evaluation of the operators Gk at the solution m] 

9. end for 

10. for all k in {l,...,d} do 
Compute and store a& (ui,ui) 
Compute and store (Gfciti, (too)v 
for all Z in {1, d} do 

Compute and store (GkU±, Giu±) v 
end for 
end for 

17. while max (£2(1^)) > tol do 

is. Find Mjv+i = ar S max (£2 (/■*)) [Selection of the parameter that maximizes the error estimator] 



[First coefficient of the reduced matrices Ak] 

[Linear terms of £ 2 ] 

[Quadratic terms of £2] 



Set ^select — ^select U {^jv +1 } 

Comp' 
Set Vi 



[Update Of Select] 



Compute Ujy +1 as the solution of E^^ 



JV+l 



[Enrichment of the reduced basis] 
[Update of the reduced right-hand side B] 



Compute and store b (^u^ +l ^j 
for all k in {1, d} do 

Compute G)~u^ +1 as the solution of Eg^ +1 k [Evaluation of the operators Gk at the 
solution Ujy +1 ] 
end for 

for all % in {1,...,N + 1} do 
for all k in {1, d} do 

Compute and store ak (ui,u^ + ^j and (uft +1 ,Uj?J 

matrices A/-] 

Compute and store (Gfctij, Goo)v 
for all I in {1, d} do 

Compute and store ^GkUi, GiUfj + ^j 
end for 
end for 
end for 

N <— N + 1 [Update the size of the reduced basis] 

end while 



[Update of the reduced 
[Linear terms of £2] 
[Quadratic terms of £2] 



Definition 2.1.1. The relative error of the computation x ® y is defined, for x (g> y 7^ 0, by 

(x®y)-(fl(fl(x)®fl(y))) 



x <8> y 



(16) 
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and the absolute error of the computation x ®y is defined, for all x®y, by 

\(x®y)-(fl(fl(x)®fl(y)))\. (17) 

If x ® y 7^ 0, a loss a significance occurs when the relative error of the computation is much 
larger than the machine error e, and ifx®y = 0, a loss of significance occurs when the absolute 
error of the computation is much larger than the machine error e. 

A well-known case of loss of significance is a sum of terms where the result is much 
smaller than the intermediate computed values. More generally, a naive implementation of 
an algorithm may lead to such a loss of significance. In such a case, one should modify the 
algorithm to ensure that each step avoids this phenomenon. In some cases, simply changing the 
order of the operations improves the situation. As an illustration, consider x = 1, y = 1 + 10 -7 
and the operation x 2 — 2xy + y 2 . This is a sum of terms where the first intermediate result in 
the sum is 14 orders larger than the result. Therefore, a loss of significance is expected. The 
relative error of this computation is about 8 x 10 -4 3> e. Computing (x — y) 2 , which is the 
factorization of the considered operation, leads to a relative error of about 10" 9 . Thus, the 
terms of the sum are only 7 orders larger than the results, leading to a less catastrophic loss 
of significance. The remedy consists in carrying out the sum before the multiplication. 

In the reduced basis context, the evaluation of the formula £2 suffers from such a loss of 
significance, as we now explain. The remedy follows the same idea as the reordering of the 
operations in the above example. 



2.2 Validity of the formulae S\ and 82 for computing the error estimator 



Consider the two formulae (12) for £\ and (14) for £2, for computing the error estimator 



Definition 2.2.1. The formula £\. is said to be valid for computing the error estimator in 
Algorithm^ with tolerance tol if 

max (£ fc (/i)) < tol. (18) 

M^^sclcct 

From a theoretical viewpoint, the residual G^u^ vanishes V/U E Select- Hence, any formula 
for computing the residual-based error estimator vanishes as well and therefore is valid with 
any tolerance. However, the validity of a formula for computing the error estimator is to be 
considered in the presence of some adverse phenomenon introducing errors in the computation, 
such as round-off errors and inexact reduced basis functions. In this context, the criterion 
max > tol clearly indicates that this adverse phenomenon prevents the convergence 

"^select 

of Algorithm [T] since the stopping criterion cannot be reached by increasing ^select j see Figure 

m 

We examine the validity of the formulae £\ and £2 for computing the error estimator in 
the presence of two independent phenomena: round-off errors and approximate reduced basis 
functions u% (in the context of inexact linear algebra solvers for E^). 



2.2.1 Round-off errors 

In what follows, we do not try to quantify the accuracy of the solution to E^, but how accurate 
the two formulae £ 1 and £2 are for computing the error estimator. We thus neglect the round- 
off errors introduced when solving E^ and E^, so that the vectors of the reduced basis Uj and 
the reduced solutions are considered free of round-off errors. 
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parameter value 




parameter value 




Figure 1: Schematic illustration of Definition 2.2.1 with Select = /X4}. Left: the 



formula £ k is valid for computing the error estimator in Algorithm [T] with tolerance tol ; right: 
the formula is not valid. 



Proposition 2.2.1. Let p G Select and maxfl(£k(/i)) denote the evaluation of the estimator 
using the formula £ k , when the maximum accumulation of round-off errors occurs. There 
holds 

maxfl(£L( M )) > 2p~ 1 6e, 
maxfl(£ 2 (//))> 2/3- 1 <5 v / e, 
where j3 m = sup (/3^) and e is the machine precision. 

Proof. Let \i G V se \ect- First, suppose that YliLi Efc=i lii^k^GkUi and 2Re(s t x^)+x*^Sx l _ l 

are computed free of round-off errors. We then have Ei=i Ylt=i li(l J ') a k{p)GkUi = —Goo 
and 2Re(s*x At ) + xfiSxp = —5 2 . We denote ((fi)i<i<N the basis functions of the Galerkin 
discretization of En. For all 1 < p < N, we obtain 



N d 

\fl((Goo) P ) + f ji(/j)a k (fj,)(G k Uijp 1 



i=l k=l 

N d N d 

< \ fl((G o)p) ~ (Goo)pl + \fl(^2^2ji(fi)ak(fi)(GkUi)p) -^^1i(( i ) a k(t i )( G kUi)( 

i=l k=l i=l k=l 

< 2(G o) P e. 

^^v 



Therefore, maxfl(/Z((G o)p)+/KEiIi ELl 7^K(/u)(G fc ^) p )) = 2(G 00 ) P e. Likewise, maxn(/7(5 2 )+ 

TV 



fl(2Re(s t x IM ) + x^Sx^)) = 25 2 e. Now, we consider that the intermediate summations in 



Ej=i Efc=i 7i(l i ) a k(^)GkUi also produce round-off errors ; we end up with a lower bound for 
maxfl(£i(/i)) and the result follows. We conclude likewise for £2. □ 



Remark 2.2.1 (Round-off errors). In "practice, in the computation of the sums Ei=i Efc=i li{l 1 ) a k{l 1 )GkU 
and 2Re(s t x fl ) + x**/5x M; £/te magnitude of the largest intermediate result is respectively of the 
order of Goo and 5 . We indeed observe in our simulations that the round-off errors on £\ scale 



2.2.1 



like e, while the round-off errors on £2 scale like y/e (see Section 3.4). From Definition 
if we suppose that, V/z G V se \ect, £\{p) < 2/3 rn 1 5e and £2^) < 2/3~ 1 5e 2 , then the formulae 
£\ and £2 are valid for computing the error estimator in Algorithm [7] with tolerance tol if, 
respectively, 

for£i, 2Bl. 1 Se < tol, 

7™ ~ (20) 

for£ 2 , 2/?- 1 5 v / e < tol. 
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As explained at the end of Section |2.1[ the computation of a polynomial using a factorized 
form is more accurate than using the developed form, in particular at points close to its roots. 

Here, (^£2^)) i s a multivariate polynomial of order 2 in x M computed in a developed form, 
whereas the scalar product (G^u^, G^u^y is not developed in the formula £\. The developed 
formula £2 was introduced to get an efficient online formula. We will see in Section [3] how 
efficiency and accuracy can be both preserved using another rewriting of the error estimator. 

Remark 2.2.2 (Accuracy levels). One might object that accuracy levels such as fi^Se can 
only be reached in simple academic problems, whereas in more complex problems, Algorithm^ 
is performed with a larger value for tol. However, as underlined in I33\j . accuracy issues can 
happen in non-compliant and nonlinear cases. Accuracy is also important in RB procedures for 
real-time control, where the feedback loop depends on an efficient error estimator. Finally, in 
industrial studies, computations of large problems are often carried-out in simple precision to 
save some computational time. In this context, it is preferable that the online error estimator 
is computed with the same level of accuracy as the solutions to E„. 

Remark 2.2.3 (Improved floating point arithmetic). Increasing the machine precision from e 



to e (quadruple precision) for computing the coefficients in (14), as well as for the evaluation 
of the multivariate polynomial in x^, is a first solution. There are also methods allowing one 
to double the precision of the evaluation of a polynomial while keeping the double precision for- 
mat, namely compensated schemes. For instance, the compensated Horner scheme in double 
precision \21^ doubles the precision and is faster than the full quadruple precision implemen- 
tation. However, this corresponds to representing the result of the intermediate operations by 
two doubles, one for the value in double precision and another one for the subsequent digits. 
Therefore, these strategies are equivalent to quadruple precision ( except for the computational 
savings in evaluating the error estimator). Moreover, since current architectures are optimized 
for the double precision format, changing the floating point arithmetic can potentially degrade 
software performance. 

Remark 2.2.4 (Goal-oriented case, round-off errors). The same analysis can be carried-out 
in the goal-oriented case. Let /j, £ Select- There holds 

maxfl(f 1 so ( / x))>2fe)" 1 <5 7 e 2 , 

_! (21) 
maxfl(ff Gu)) > 2 (fa) Sje, 

where fa = inf (j3~) and 7 := [|Q||y- We observe as well in our simulations that the 

round-off errors on £f° scale like e 2 , while the round-off errors on £f° scale like e (see Section 

4). If we suppose that, V/z G V se \ect, £i°(aO — 2 (Pm) &l e2 an d ^Til 1 ) — 2 ($m) ^7 e > ^ e 
formulae £f° and £f are valid for computing the error estimator in AlgorithmUj with tolerance 
tol if, respectively, 



forff 5 , 2 (fa) ^e 2 <tol, 



forff, 2(M) 5 7 e < tol. 



-1 



(22) 



2.2.2 Approximate reduced basis functions 

In practice, in particular for large-scale simulations, the accuracy of the RB procedure is also 
limited by the numerical method used for computing the RB functions. We recall that for a 
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given value /i G Select! En consists in solving a linear system of size N of the form A^IJ^ = B. 



Thus, for n G Atrial, the formulae £1 (12) and £2 (14) for the error estimator are based on the 



computation of the residual of for the reduced solution u^: ||G^u^||v = ||A^f7^ — -B||*y, 
where V $ G M. N , ||$||*v = sup ^'^"^ , where we recall that ((fi)i<i<N is the basis 

functions of the Galerkin discretization of E^, see |1 T|, §9.1.5]. In what follows, we suppose 
that the formulae £\ and £2 are free of round-off errors (therefore, V/i G Atrial , £i(a*) = £2 (/•*)), 
but the problem E 1 ^ is not exactly solved, leading to approximate reduced basis functions such 
that the residuals do not vanish. Hence, for all \i G "P S eiect> £l(/-0 = £2^) > 0. 

Proposition 2.2.2 (Approximate reduced basis functions). If the reduced basis functions are 
computed using an iterative solver with the following stopping criterion on the normalized 
residual 

^ € Atrial, \\ A ^ -B\Uv> < ^ (23) 

||±> ||*V' 

then the formulae £\ and £2 are valid for computing the error estimator in Algorithm^ with 
tolerance tol if 

fatSZ < tol. (24) 



Proof. Let k G {1,2}, let /x G "Pseiect and suppose that the stopping criterion (23) is satisfied. 
Then, = u^, but u M does not exactly solve E^. First, by definition of the || • ||*v norm, 

\\B\Uv = sup *Sn^ < \\b\\v = HGoollv = Then, ||G M fi M [| v = sup (G ^jf ;" )v = 



sup ^XZ = SU P v 7 = \ ~ B W*y- Therefore, 



Nlv Vf £ N HEfciVwIlv 
£ fc (/i) = ^HG^IIv = /S-^I^C^-BIUv' = ^WA^-BIU < p- x \\B\\^ < < fa 

Hence, if (3^5^ < tol, the result follows from Definition 



2.2.1 



□ 



Since the || • ||*y norm is hard to compute, the stopping criterion (23) is not the one used 
in practice in linear algebra solvers. 

Remark 2.2.5 (Goal-oriented case, approximate reduced basis functions). The formulae 
£f° and £f° are valid for computing the error estimator in Algorithm [7] with tolerance tol 

*/(^)~W 2 <tol. 
2.2.3 Synthesis 

Taking into account the round-off errors in the computation of the error estimator and the 



stopping criterion of an iterative solver, and supposing that the bounds (19) and (21) are 
reached, the formulae £\ and £2 are valid for computing the error estimator in Algorithm [T] 
with tolerance tol if, respectively, 

for £1, 2/3 m 1 <5 max (f , e) < tol, 

for £2, 2/3~ 1 <5max (£, y/e) < tol, 

and the formulae £f° and £ f° are valid for computing the error estimator in Algorithm [l] with 
tolerance tol if, respectively, 



for £f°, 2 ( Pi ) 5j max (£ 2 , e 2 ) < tol, 

V _! (26) 
for £f°, 2 5-y max (£ 2 , e) < tol. 
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3 New procedures for accurate and efficient evaluation of the 
error estimator 



The formula E\ for computing the error estimator is accurate but non efficient, whereas the 
formula £2 is efficient but non accurate. In this section, efficient methods achieving the same 
accuracy levels as £\ are devised. 



3.1 Procedure 1: rewriting £ 2 

We first present the procedure proposed in \9\. We consider that a reduced basis of size N 
has been constructed. Let a := 1 + 2dN + (dN) 2 . For a given [i G Atrial an d the resulting 
Up G Spanjui, u^} solving Q, we define X(n) G C a as the vector with components 
(1,^,^^), where x^j = a>k{[i)"li{pi) (we recall that 7«(/i) are the coefficients of the 
reduced solution in the reduced basis and «&(//) the coefficients of the affine decomposition of 



a M in (10)), with 1 < I, J < dN (with I = i + N(k - 1), such that 1 < i < N, 1 < k < d, and 



with J = j + N(l — 1), such that 1 < j < N, 1 < I < d). We can write the right-hand side of 



(14) as a linear form in X{ix) as follows 



5 2 + 2Re{s t x M ) + xpx^ = ^ ^X p (/i), (27) 

P =i 

where t p is independent of \i and X p (fi) is the p-th component of X(fj,). 

Now, in the offline stage, we take a values (think of random values) /x r G "Ptriai, t G 
{1, ...,<r}, of the parameter /x. Then, we compute the vectors X(fi r ) and the quantities 

V r :=Y,tp x P^r)- (28) 
P =i 

Finally, we define T G C CT><<T as the matrix whose columns are formed by the vectors X(/j, r ), 
that is, Tpr = Xpd^r) for 1 < p, r < a. We assume that T is invertible, which always happens 
to be the case in our simulations. 

Then, in the online stage, suppose that we want to evaluate the error estimator for the 
RB solution computed at the parameter value fi G "Ptriai- We compute the vector X(fi) 
and solve the linear system 

T\(j i )=X( Ji ), (29) 
for A(/z) G C 7 . We then obtain X(fi) = Yf r =i K(^)X(fi r ) and 

a cr a 

y]t p X p (fi) = t p \ r {n)X v (n r ) = ^A r (//)V;. (30) 

p=l p,r=l r=l 

This yields the new formula for computing the error estimator, 

^-^(tM^r) • (31) 



\r=l 



The quantities V r are precomputed as V r = (^(3 /Ir £i(fj ir )j . Computing £3 requires solving 

TX(fi) = X(fi) and summing the a precomputed quantities. Since the complexity of this 
procedure is independent of N, the error estimator £3 is online efficient. 
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Notice that £i(/u), £2(1^), and £3^) are equal in exact arithmetic. Following the idea of 
reordering the operations of a formula to modify its behavior in regard of round-off errors, 
the formula £3 involves the evaluation of a multivariate polynomial of order two as a linear 
combination of accurate values of this polynomial at different points, whereas £2 consists in 
evaluating this polynomial from its coefficients. In the formula £2, the first term of the sum, 



8 , is fixed. Hence, as we saw in Remark 2.2.1 computed values for £2 (/•*)> M £ ^select; depend 



on 8. In £3, since the quantities V r correspond to values of the error estimator for parameters 
in VtriaXi the terms in the sum of (31) can be of low magnitude. More precisely, if the size of 
the reduced basis is large enough, the quantities V r can be much lower than 8 2 . In such a case, 
the maximum intermediate result in the sum of formula £3 is much lower than the maximum 
intermediate result in the sum of formula £2 (which is lower bounded by 8 2 ). Therefore, the 
formula £3 can be valid for lower tolerances than the formula £2. 

As pointed out in [9], the matrix T exhibits in practice large condition numbers, and there 
is no guarantee that it is actually invertible. We will see in Section [4] for a three-dimensional 
acoustic scattering application that £3 can be in practice as ill-behaved as £2. Moreover, there 
is no a priori method for selecting the parameter fj, r for which the quantities V r have been 
precomputed. In the next section, we propose a new procedure that solves these problems. 



3.2 Procedure 2: improvement on Procedure 1 using the empirical inter- 
polation method 

In the formula £3, a potentially ill-conditioned problem TA(/i) = X(fi) is solved in order to 
exactly represent X(fi) by the linear combination Ylr=i ^r(n)X{n r ). Following a suggestion 
by Patera [32], we propose to approximate X(fj,) by means of an interpolation procedure. We 
want to modify the formula £3 by an interpolation formula relying on a better conditioned 
linear system. The price to pay is that the new formula £4 will not equal £\ in exact arithmetic. 
We observe in our numerical experiments that the error introduced by this interpolation is 
lower than the ones due to the poor conditioning of T in the evaluation of £3. We also look for 
a way to choose the parameters [i r for which the quantities V r have to be precomputed. We 
refer to these values for \x r as "interpolation points", and to the set of these points as "Pinter- 
Consider the function of two variables (p, /i) — > X p (fi), for p G {l,...,cr} and fi G Atrial- 
We look for an approximation of this function in the form 

V/i G 7W,Vp G {1,...,ct}, X p (h) « ^Af(M)^(M, (32) 

r=l 

where o < a. The empirical interpolation method (EIM) (more precisely the discrete EIM 
since p is a discrete variable) provides a numerical procedure to construct this approximation 
and to choose the interpolation points (see [3} l2"9"]). 

For completeness, we briefly recall the EIM and adapt the notation of [29] to the present 
context. The EIM is an offline-online procedure. During the offline stage, a basis functions 
are computed, denoted qj : Atrial 9 ^ <&'(/■*) S C, 1 < j < a. These basis functions will 
be used in the online stage to carry-out the interpolation. We define q a as the vector-valued 
map Atrial 9 /J !-> (f(^) '■= (<7i(M))i<7<CT 6 During the offline stage, a interpolation points 
fjL r G Atrial are also selected; these points are collected in the set Tenter- During the online 
stage, the matrix B a G <£ a,a , where Bfj = qi(nj), for 1 < i, j < a, is constructed. Letting 
H G "Ptriah w e solve for X a (fJ-) G C CT such that 

B*A*fa) = <?(jj), (33) 
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and compute the rank-o" interpolation defined as follows. 

Definition 3.2.1. Let 1 < k < a, the rank-k interpolation operator I k is defined such that 

k 



I k X{») :=^A r fe (/x)X(/x r ), 



(34) 



r=l 



where X k (fi) G C k solves B k \ k {^) = q k (n) 



The formula X p (fi) (I a X) p (fj,), V/z G PtriabVp G {!>•••> o"}> provides the interpolation 
formula searched for in ( 32 ) . 



Definition 3.2.2. The residual operator 5 a is defined by 

(T := Id - f. 



(35) 



Algorithm [2] presents the construction of the function q a by a greedy algorithm during the 
offline stage. 

Algorithm 2 Offline stage of the EIM 



1. Choose a > 1 

2. Set := 1 



3. Compute pi := argmax [|-Xp(')[l*~(P tri ai) 

pe{i,...,s} 

4. Compute //1 := argmax|X pi and set Pinter = {a*i} 



5. Set qi(-) :- 



Set 5 



-Xpi(-) 
Xpi (pi) 
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1 



7 while k < a do 

«. Compute pfc+i := argmax || (<5 A: ^)p 1 (-) IU°°(7^ trial ) 
pe{i,...,s} 

9. Compute Hk+i '■= ar g max l(^ fc ^)p fc+ i(A t )| 

10. Set Pinter := Pinter U {^fc+l} 

11. Set %+i(-) := 



112. 

13. 



B 



fc+i ._ 



(<5 fe ^WiOfc+i) 



k <r- k + 1 
end while 



[Number of interpolation points] 



[First interpolation point] 

[First basis function] 
[First B matrix] 



[(k + l)-th interpolation point] 
[Update of Pinter] 
[(k + l)-th basis function] 

[(k + l)-th B matrix] 
[Increment the size of the interpolation] 



Definition 3.2.3. The new formula for computing the error estimator is 



(36) 



v r=l 



where A CT (^) is the solution of (33). We recall that V r = \J3^ r Ei{^i r ) 
Proposition 3.2.1. The computation of the formula £4 is well defined and online efficient. 
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Proof. Owing to [29, Theorem 1], the matrix B is lower triangular with diagonal unity. Hence, 
det B = 1 and B is guaranteed to be invertible. The online procedure of EIM, consisting in 
solving a linear system defined by the matrix B, is thus well defined. Then, since the EIM 
procedure in carried out on X p (/i), \/p G {1, ...,cr}, V/i G "Ptriah all the computations involved 
are of complexity independent of N, even the offline part of the EIM. Finally, the complexity 
of the online part of EIM only depends on a. □ 

Remark 3.2.1 (Stopping criterion in Algorithm [2]) . For ease of presentation, we chose a 
simple stopping criterion based on a a priori fixed maximum number of interpolation points. 

Remark 3.2.2 (Complexity of the offline stage in EIM). If the cardinality of Vtxi&i is large, 
the complexity of the offline stage of the EIM can be high. 



3.3 Procedure 3: improvement of Procedure 2 using a stabilized EIM 

In practice, round-off errors are accumulated during the loop in Algorithm [2j and when search- 
ing for an accurate interpolation of X(n), the matrix B becomes non invertible at some point, 
if we keep increasing the number of interpolation points. The coefficients of the matrix B 
suffer from round-off errors, so that the relation det(B) = 1 no longer holds. To solve this 
problem, we now propose a numerical stabilization of EIM based on the following property: 

Property 3.3.1. There holds 

Vi < j, P o F = P, (37) 



where the interpolation operators P are defined by (34). 

Proof. Using [29], Lemma 1], PX G Span (gi, qj) and Pv = v for all v G Span (q%, qj). 
Therefore, P o PX = PX. □ 

In our numerical experiments, we observed that, as the number of iterations of the greedy 



procedure for the EIM grows, the relation (37) is no longer verified numerically, due to accu- 



mulations of round-off errors. These numerical instabilities can be compensated in the same 
fashion as the Gram-Schmidt orthonormalization process is stabilized in practice (see [211 
chapter 5.2.8]). The Gram-Schmidt algorithm transforms a linearly independent family of 
vectors {vj} into an orthonormal basis {u{\. To simplify the presentation, we suppose in what 
follows that the normalization step is not carried out. Consider the orthogonalization step for 
the /c-th vector. We denote Ii k the projection operator on Span(ui, u^), and 5 k := Id — II fc . 
For EIM, we suppose that (k — 1) interpolation operators P, 1 < i < k — 1, have been con- 
structed, and we wish to construct the k-th interpolation operator I k . A comparison between 
the stabilized Gram-Schmidt orthonormalization procedure and the proposed stabilization for 
EIM is presented in Table [T} 

Proposition 3.3.1. Let k G N*. In exact arithmetic, the following relations hold for the 
residuals defined in Table^^ 5 k tab v = 8 k v. 



Proof. We propose to prove by recursion that, Vi < k, = 5 l . The case i = 1 is clear 
from the definition of the first intermediate residual in Table [TJ Let i < k and suppose that 
S^- 1 = Id - P- 1 for EIM. There holds 

^: b = S 1 -^ S 1 =Id-/ i - 1 -^ + / i o^ 1 = Id-/ l = ^, (38) 



since P o P^ 1 = P~ L owing to Property 3.3.1. The results follow from the case i = k. The 
same relation is proved likewise for Gram-Schmidt, for which II* o IF -1 = IF -1 holds as well 
by definition of the projection operators II 1 . □ 
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stabilized Gram-Schmidt 


stabilized EIM 


global input 


(«Xj -;Vs-) basis of C CT 


V : Atrial -> 


classical residual at step k 


S k v k = v k - U k v k 


(5*«)(/ 1 )=«( M )-(J*t;)(/i) 


intermediate residuals at step k 


s stL v k = v k - nV. 

ck.2 jfe,l tt2x*:-1 

ck,k ck.k— 1 Tri; irk,k— 1 
<> B tab t >* = <> S tab U fe - 11 <> S tab V k 


(^)(a*) = «(m)-(/ 1 «)(m) 

0&ab")(M) = (^>)(i")-/ 2 (4'ab«)(M),--. 
(^b«)(M) = (&M(m) ~ ^aTM(M) 


stabilized residual at step k 


%ab U fc = ^tab^fe 


= (&)(/*) 


global output 


(^stab^l^stab^.-^ftab^) 

orthogonal basis of Span(i>i, v&) 


(J*«)( M ) 



Table 1: Comparison between stabilized Gram-Schmidt and stabilized EIM. 



Definition 3.3.1 (Stabilized E IM) . The stabilized EIM consists in the same offline procedure 
as the one described in Section 3.2, except that the residuals 5 k are replaced by the stabilized 



residuals <5g tab defined in Table |i| The online stage is the same as the classical EIM. 

The stabilized Gram-Schmidt generates a set of vectors much less polluted by round-off 
errors (see jHEi])- By analogy we expect that the stabilized EIM produces a more accurate 
interpolation procedure than the classical EIM that is, much less polluted by round-off errors. 
This is numerically verified in Figure [2J where det(-B CT ) and cond(i3 <T ) are represented as a 
function of a. If the method is stable, then det(B a ) = 1 should hold throughout the process. 
Anticipating on the numerical applications, we consider the test case described in Section 
namely a one-dimensional diffusion problem for which N = 7, d = 2, and a = 225. In 



3.4 



Figure[2jis plotted det(-B CT ) as a functions of a. The condition number is also plotted to see how 
numerical instabilities worsen the conditioning of the matrix B a . The stabilized EIM behaves 
as intended. The classical EIM curve stops since the matrix B a becomes non invertible at some 
point: a parameter already in Pinter has been selected by the greedy algorithm. Invertibility 
can artificially be recovered by ensuring that the new interpolation point cannot be an element 
of the current set Pinter- We call this procedure EIM with unique choice. However, this fix is 
not completely satisfactory, since det(B a ) = 1 is not satisfied. Moreover, cond(-B CT ) is more 
ill-behaved with this procedure than with the stabilized EIM. 

Remark 3.3.1 (Computational cost and variant of stabilized EIM). The computational cost 
of the stabilized EIM is much more than that of the classical EIM, since the stabilized residual 
requires as many calls to a classical residual as the number of selected interpolation points. 
One can think of a cheaper procedure by monitoring det(B a ) and adding some intermediate 
residuals S^A, until det(B a ) is close enough to 1. 



3.4 Summary and illustration 

The advantages and drawbacks of the four considered formulae for computing the error esti- 
mator are summarized in Table [2j 

To illustrate, consider as in [9] a one-dimensional linear diffusion problem, namely the 
boundary value problem — u" + fj,u = 1 on ]0, 1[ with u(0) = u(l) = 0, where \i G V := [1, 100] 
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Figure 2: Determinant (left) and condition number (right) of the matrix B a as a function of 
a, for classical EIM, classical EIM with unique choice and stabilized EIM. The classical EIM 
curves stop at 21 interpolation points since B a becomes non invertible at 22 points. 



Property 


Si 


£2 


£3 


£4 


Online efficient 


No 


Yes 


Yes 


Yes 


Unconditionally well-posed 


Yes 


Yes 


No 


Yes 


Dependence on e of the observed accuracy 


e 


V~e 


e, if well-posed 


e 


Equals £\ in exact arithmetics 




Yes 


Yes 


Yes, if a = a 
No, if a < a 



Table 2: Comparison of the considered formulae for computing the error estimator. 



is the parameter. The analytic solution is 



u(x) = (cosh {JJix) — 1) -| (v^f) — siuh (J~jlx) . (39) 

/U //sinh(- v //^J 



The Lax-Milgram theory is valid, and the coercivity constant is bounded below by 1 in the H 1 - 
norm. The error estimator is given by £i(/u) = HG^M^H^-ino^). Lagrange Pi finite elements are 
used with uniform mesh cells of length 0.005. The set 'Ptrial consists of 1000 points uniformly 
distributed in [1, 100]. The RB method is carried out until the formula £2 suffers from round- 
off errors, which already happens for a reduced basis of size N = 7 (since d = 2, we obtain 
a = 225). A direct solver is used, so that we only expect round-off errors in the formulae for 
computing the error estimator. 

In Figure [3j we see that the classical formula £2 is not valid for computing the error 
estimator in Algorithm [l] with any tolerance below 10~~ 7 , whereas the formulae £\, £3 and £4 
are valid with tolerance down to 10 -14 . The difference is of 7 orders of magnitude ; given that 



y/e ~ 10 ; this is consistent with Remark 2.2.1 and Section 3.1 



In Figure |4j we observe that instabilities occur in the formula £3, especially for parameter 
values close to the elements of Select- This is due to the poor conditioning of the matrix T 



17 




when solving (29). Therefore, the new formula £4 based on the EIM introduces less numerical 
errors than £3. 



7 
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Figure 4: Comparison of the error estimators £3 and £4, with respect to the reference £\. 
In Figure |EJ we plot niax (£4(^4)) as a function of <r. From this figure and Definition 



2.2.1 



select 



we deduce that for a < 23, the choice of any tolerance larger than 10 -12 will make the 
formula £4 valid. If we want to consider a tolerance of the order of 10 -14 , we need a > 23. In 
this example, the decrease of max (£4(^1)) with respect to a is very steep at a = 23. We will 

^select 



see in Section |4 an example were max (£4(^1)) decreases more slowly. Namely, on Figure 10 

/iS'Psclcct 

from left to right, max (£4(11)) equals successively 5 x 1CP 8 , 1CP 8 , 8 x 1CP 9 and 2 x 10~ 9 , 
for a = 20, 30, 40 and 50.°* 
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Figure 5: max (S^p)) as a function of a. 

oct 

4 Application to a three-dimensional acoustic scattering prob- 
lem 

4.1 Formulation of the problem 

We consider a ball Q l C M 3 with boundary T and Q e := M 3 \0 J , see Figure [6j We consider 
a monopole source located in fi e . The surface of the ball is impedant, meaning that any 
incident wave will be partially absorbed and partially scattered. The proportion of absorbed 
and scattered parts are quantified by the impedance coefficient fj,, which is used in a Robin 
boundary condition at T. We are interested in the computation of the scattered field p sc in Q, e . 
We denote p mc the known pressure field created by the source in the absence of the sphere; 
the total acoustic field in f2 e is the sum of p mc and p sc . 

acoustic monopole 




Figure 6: Geometry for the three-dimensional acoustic scattering problem 

We define the distribution v : Q e U — > C such that v\ ni = —pine, ^|n e = Psc- We denote 
A and \ the jumps of the Neumann and Dirichlet traces of v across T. The Robin boundary 
condition writes A + — x = 0. Since v solves the homogeneous Helmholtz equation in fl e and 



19 



in Q l and satisfies the Sommerfeld radiation at infinity, there holds 



V = -S\ + V X inft e Uir, (40) 
where S and T> are respectively the single- and double-layer potentials. Taking the inte- 



rior Dirichlet and Neumann traces of v in equation (40) and injecting the Robin boundary 
condition, we obtain 



N 



ik j 
2/i J 



D 



D 



-S 



IE T 
2k 1 













X 




7l Pine 




A 




_ -7o>inc _ 



(41) 



where k is the wave number of the monopole source, N, D, D and S are classical boundary 
integral operators (see [US]), and 7(7 Pine and 7{~ p mc are respectively the interior Dirichlet and 
Neumann traces of the known function Pi nc . Solving one of these two equations, together 
with the Robin boundary condition, is sufficient. The software we are using, ACTIPOLE (see 



161 [38] ) , deals with the block system defined in (41), which presents the advantage of being 



invertible for all frequencies of the source, even when the surface is only supposed Lipschitz. 



We denote the block operator defined by the left-hand side of (41). From [251 EH 136] . 
we infer that A^ is a bounded bijective operator from H^(T) x L 2 (F) into H~2(T) x L 2 (F) 
(see also [10 ). The variational form is as follows: find (x, A) G H2(T) x L 2 (T) such that 
V(x,A) Gj Ffl(r)xL 2 (r), 



t / ik 



A,SA + |a 



A, 70Pir 



(42) 



denotes the Hz(T) x H 2 duality product and < 



> denotes the L 2 (T) inner 



where (•, 
product. 

Let and be respectively the Lagrange finite element spaces Pi and Fq on T. Let 
(<j>i)i<i<p and (ipj)i<j<P' be the usual bases of and V h ° of size P and P', respectively. 
The product space x V® is a conforming approximation of i?2(r) x L 2 (r). The numerical 
resolution is carried out with a Galerkin procedure on x using the boundary element 
method (BEM). From |25j . the obtained discrete approximation of the problem (42) is inf-sup 
stable, as required by the RB method (see also [TO]). 



4.2 Application of the RB method 

We take as parameter for the RB method the value of the impedance which is supposed 
here to be a positive real number. To recover an affine dependence on the parameter fj,, we 
write the BEM matrix in the form A^ = ai{^)A\ + 02(^)^2 + 03(^)^3, so that d = 3 in the 



affine decomposition (10) with aiQu) = 1, «2(^) = 77 and as(fi) = fi. Specifically, 







1 < i < P 




1 < i 


< P 


A x = 




1 < j < p 




l<3 


< P' 




1 < i < P> 
1<3<P 




1 < i 
1 < j 


< P' 

< P' 



(43) 
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-f (<l>i,<f>j) 


1 < i 


< P 


(0) 


1 < i 


< P 


A 2 = 




1 <3 


< P 




1 < 3 


< P' 


(0) 

i < 


i < P' 




(0) 


1 < i 


< P' 




1 < 


3 < P 






1 < 3 


< P' 



(0) 



l<i< P 
1<3<P 



(0) 



1 < i < P' 
1<3<P' 



(0) 



1 < i< P 
1<3<P' 



4 (V'i,^ 



i < i < p' 

±<3<P' 

(44) 

In the general-purpose RB, the quantity of interest is the pair of potentials (x, A) on T. 
For the goal-oriented case, we consider the value of the pressure at a given point of space. If 
this point if far enough from the sphere, approximations can be made in the representation 
formula for the pressure. This is the far-field approximation, which consists in a linear form 
Q acting on the solution pair (x, A) as Q ((x, A)) € C 2 , with 



Q(x,A) 



V 



-ik- 



-ife||a;|| 



ik 



4:ir\\x\\2 

e -ik\\x\\ 2 



-iky- 



n(y),x(y) 



4-7T \\x\ 



-iky- 



) 



6 



(45) 



For simplicity of the computation, we will take the Euclidian norm on the vectors of C p+P ' 
instead of the H2(T) x L 2 (T) norms of the reconstructed functions. This way, the Riesz 
isomorphism J is simply the identity. Therefore, the computation of the terms G^u^, as well 
as that of the terms G^Ui, does not require to invert the stiffness matrix as in (15). 

We define two test-cases: (i) one impedant sphere (d = 3), with N 
[0.9, 1.1], (ii) two impedant spheres (d = 5), with N = 1561 and /j, £ V :- 



584 and \i G V : = 
[0.99, 1.01] 2 . 

We present visualizations of the diffracted pressure field, at a random parameter value ji, 
for the test-cases (i) with #7-triai = 100 and N = 10 in Figure [TJ (ii) with #'Ptriai = 225 and 
N = 10 in Figure M 
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Figure 7: Real part of the pressure field for a) the BEM solution, b) the RB solution (basis 
of size 10) and c) the difference 



4.3 Error estimator curves 

We present the error estimator curves for the test-cases (i) general-purpose RB, ^'Ptriai = 100, 
(N,a,a) = (2, 7, 49), (3, 10, 100), (4, 20, 169), and (5,30,256), in Figure R ; (ii) goal-oriented 
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Figure 8: Real part of the pressure field for a) the BEM solution, b) the RB solution (basis 
of size 10) and c) the difference 



RB, #P t riai = 225, N = 8, a = 60, and a = 1681, in Figure 



11 



In the experiment (i), we investigate the behavior of £4 when a increases in Figure 10 



Although the four cases lead to the same local maxima, we see that we need more interpolation 
points a as we want the formula £4 to be valid for smaller tolerances. 

In the experiment (ii), the curve for £3 is quite poor, and we do not observe the level 
of accuracy we observed in the previous simple examples. Here, the matrix T (29) is so ill- 
conditioned that the numerical errors introduced by its resolution are larger than the ones 
introduced by the formula £2- The classical formula £2, in this configuration, requires 1681 
offline precomputations, which is much larger than the 60 quantities V r used to compute 

we see that argmax (£4^)) = (1,1) and £"4(1,1) ~ 10~ 16 ; 

^select 

60 is valid for computing the error estimator in Algorithm 



the formula £4. From Figure 



11 



therefore, the formula £4 with a 
1 with tol = 10~ 16 . 



Conclusion 

In this work, we have extended the ideas of [9] by proposing a more stable numerical procedure, 
using the empirical interpolation method, to represent the a posteriori error estimator in 
the reduced basis as a linear combination of its values at given parameter values, called 
interpolation points. Moreover, this method provides a way of choosing the interpolation 
points, and yields better accuracy levels than the classical a posteriori error estimator and than 
the procedure proposed in [SJ. Besides, our new procedure may require less precomputations 
than the classical a posteriori error estimator. The new error estimator derived herein can be 
of particular interest in two situations: (i) when accurate evaluations of the error estimator are 
needed and (ii) when one wants to decrease the computational time of the offline procedure. 
However, contrary to the classical formula for computing the a posteriori error estimator, 



all the quantities V r defined in (28) have to be recomputed when adding a new vector in 



the reduced basis. Besides, the complexity of the offline stage of the EIM depends on the 
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Figure 9: Error estimator curves with respect to the impedance coefficient, with TV equal to 
2, 3, 4, and 5 (from left to right and top to bottom). The curve £2 computed in quadruple 
precision superimposes to £\. 




Figure 10: Error estimator curve £4 with respect to the impedance coefficient, with N = 5 
and a equal to 20, 30, 40, and 50 (from left to right). 



cardinality of Atrial- 
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Figure 11: Error estimator curves (logarithmic sale) as a function of the impedance coefficients: 
a) £i, b) £2, c) £3, d) £4, and e) £2 computed in quadruple precision. 
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